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Code Description 

This work uses a discontinuous-Galerkin spectral-element method (DGSEM) to solve the compressible Navier-Stokes 
equations [1—3]. The inviscid flux is computed using the approximate Riemann solver of Roe [4]. The viscous fluxes 
are computed using the second form of Bassi and Rebay (BR2) [5] in a manner consistent with the spectral-element 
approximation. The method of lines with the classical 4th-order explicit Runge-Kutta scheme is used for time 
integration. Results for polynomial orders up to p = 15 (16th order) are presented. 

The code is parallelized using the Message Passing Interface (MPI). The computations presented in this work are 
performed using the Sandy Bridge nodes of the NASA Pleiades supercomputer at NASA Ames Research Center. 
Each Sandy Bridge node consists of 2 eight-core Intel Xeon E5-2670 processors with a clock speed of 2.6Ghz and 
2GB per core memory. On a Sandy Bridge node the Tau Benchmark [6] runs in a time of 7.6s. 


Case 1.6: Vortex Transport by Uniform Flow 


Case Summary 


The isentropic vortex convection problem is initialized with a perturbation about a uniform flow given by: 
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where Uoo is the convecting speed of the vortex, /3 is the vortex strength, R the characteristic radius, and (x c , y c ) 
the vortex center. The convection velocity is angled 30° from the ^-direction so that the flow does not align with 
the mesh. The exact solution to this flow is a pure convection of the vortex at speed Uoo- 

For the unsteady simulation we employ a time-step set by a CFL condition based on the acoustic speed: At = h 
where h is the element length scale, c is the freestream speed of sound, N = p + 1 is the solution order and CFL = 0.5. 
Given this temporal resolution, the spatial error dominates, such that reduction of the time-step changes the error 
by less than 0.1%. 
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Meshes 

The code used to perform the simulation is written to be three-dimensional. As such this two-dimensional test-case 
was run with the fully three dimensional version of the code by extruding the two-dimensional mesh with a single 
element in the third coordinate direction. We report the number of degrees of freedom as the equivalent number of 
two-dimensional degrees of freedom. 


Results 

This test case has an analytical solution so that the exact error may be computed. The L 2 error was computed 
in a manner consistent with our spectral element discretization; namely, using the collocation points as quadrature 
points. The error is measured after 50 convective time periods. 


Fast Vortex, M = 0.5 

Figure 1 plots the L 2 error versus h and L 2 versus TAU work units (i.e. the CPU time normalized by the time 
required to run the Tau Benchmark [6]). Table 1 provides the same data in tabular form. As seen in Figure 1 
increasing the solution order allows for a lower error tolerance to be met with fewer degrees of freedom. In terms 
of TAU work units 4th-, 8tli- and 16th-order methods have a significant advantage over the 2nd-order method. We 
note that since we are using a three-dimensional code to perform our simulations with a single element in the third 
coordinate direction, the total number of actual (3D) degrees of freedom used in our simulation increases linearly 
with solution order. Even with this artificially inflated cost, increasing the solution order is an efficient means for 
achieving low error. 


DOF /dir 

P 

Order 

Elem/dir 

TAU 

Work 

L2 

Error 

32 

1 

2 

16 

1.37 

X 

10 1 

1.53 

X 

10~ 2 

64 

1 

2 

32 

1.03 

X 

10 2 

1.44 

X 

10~ 2 

128 

1 

2 

64 

7.98 

X 

10 2 

1.55 

X 

10" 2 

256 

1 

2 

128 

6.60 

X 

10 3 

6.49 

X 

10~ 3 

512 

1 

2 

256 

5.34 

X 

10 4 

1.65 

X 

10~ 3 

1024 

1 

2 

512 

4.35 

X 

10 5 

3.31 

X 

10” 4 

32 

3 

4 

8 

3.89 

X 

10 1 

1.75 

X 

10~ 2 

64 

3 

4 

16 

2.88 

X 

10 2 

2.16 

X 

10~ 2 

128 

3 

4 

32 

2.25 

X 

10 3 

1.14 

X 

10” 3 

256 

3 

4 

64 

1.84 

X 

10 4 

3.17 

X 

10“ 6 

512 

3 

4 

128 

1.49 

X 

10 5 

4.13 

X 

10" 8 

32 

7 

8 

4 

1.64 

X 

hF“ 

2.04 

X 

10~ 2 

64 

7 

8 

8 

1.23 

X 

10 3 

1.86 

X 

10~ 3 

128 

7 

8 

16 

9.85 

X 

10 3 

4.91 

X 

10- 7 

256 

7 

8 

32 

7.84 

X 

10 4 

6.06 

X 

10- 9 

32 

15 

16 

2 

6.26 

X 

hF“ 

2.43 

X 

10~ 2 

64 

15 

16 

4 

5.16 

X 

10 3 

5.45 

X 

10~ 4 

128 

15 

16 

8 

4.08 

X 

10 4 

1.64 

X 

10" 8 


Table 1: Case summary for the isentropic vortex convection, M = 0.5. 


Slow Vortex, M = 0.05 

Figure 2 plots the L 2 error versus h and L 2 versus TAU work units for the slow moving vortex test case. The 
corresponding tabular data is given in Table 2. As in the fast moving test case, there are clear benefits to using 
higher-order. 
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h=V DOF 



(a) L 2 Error vs h 


(b) 1,2 Error vs Work 


Figure 1: Error convergence for the isentropic vortex convection, M = 0.5. 


DOF/dir 

V 

Order 

Elem/dir 

TAU Work 

L 2 Error 

32 

1 

2 

16 

1.39 x 10 2 

1.44 x 10” 4 

64 

1 

2 

32 

1.03 x 10 3 

1.40 x 10~ 4 

128 

1 

2 

64 

8.00 x 10 3 

1.21 x 10~ 4 

256 

1 

2 

128 

6.59 x 10 4 

7.49 x 10~ 5 

32 

3 

4 

8 

3.90 x 10 2 

1.16 x 10~ 4 

64 

3 

4 

16 

2.89 x 10 3 

4.44 x 10~ 5 

128 

3 

4 

32 

2.24 x 10 4 

2.28 x 10~ 6 

256 

3 

4 

64 

2.03 x 10 5 

2.46 x KT 8 

32 

7 

8 

4 

1.64 x 10 3 

6.16 x 10~ 5 

64 

7 

8 

8 

1.23 x 10 4 

1.22 x 10~ 6 

128 

7 

8 

16 

9.84 x 10 4 

3.08 x 10" 10 

32 

15 

16 

2 

2.34 x 10 4 

1.88 x 10~ 5 

64 

15 

16 

4 

2.08 x 10 5 

1.41 x 10" 8 


Table 2: Case summary for the isentropic vortex convection, M = 0.05. 
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(a) Z/2 Error vs h (b) L 2 Error vs Work 

Figure 2: Error convergence for the isentropic vortex convection, M = 0.05. 


Case 3.5: Direct Numerical Simulation of the Taylor-Green Vortex at 

M 0 = 0.1, Re = 1600 


Case Summary 


The Taylor-Green vortex flow is simulated using the compressible Navier-Stokes equations at Mq = 0.1. The flow is 
solved on an isotropic domain which spans [0, 27rL] in each coordinate direction. The initial conditions are given by: 


u 

v 

w 

V 


Vq sin (x/L) cos (y/L) cos (z/L) 

— Vo cos (x/L) sin (y/L) cos {z/L) 

0 



1 

16 


(cos(2x) + sin(2j/)) (cos( 2 . 2 ) + 2)) 


( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 


where u,v and w are the components of the velocity in the x, y and ^-directions, p is the pressure and p is the 
density. The flow is initialized to be isothermal (p = ^ = RTq). The flow is computed at a Reynolds number of 
Re = PaX ° L = 1600, where p is the viscosity. The Prandtl number is Pr = 0.71, while the bulk viscosity is given by 
the Stokes hypothesis: A = — |/r. 

The unsteady simulation is performed for 20f c , where t c = pjj is the characteristic convective time. The time-step is 
set by the CFL condition based on the acoustic speed: At = where h is the element length scale, N is the 

solution order and CFL = 0.5. 


Meshes 

Direct numerical simulation was performed using three different mesh sizes for each polynomial order considered, 
such that the total number of degrees of freedom in each coordinate direction was 128, 192, 256 or 384. Table 3 
summarizes the mesh sizes and polynomial orders used. Table 3 also gives cost of the simulation in terms of TAU 
work units. At 8th and 16th order the DG discretization does not provide sufficient numerical dissipation on the 
under resolved cases (128 3 for both 8th and 16th order and 192 3 for 16th order) and these cases are numerically 
unstable. 
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DOF/dir 

P 

Order 

Elem/dir 

TAU Work 

128 

1 

2 

96 

3.30 x 10 4 

128 

3 

4 

48 

4.56 x 10 4 

128 

7 

8 

24 

unstable 

128 

15 

16 

12 

unstable 

192 

1 

2 

96 

1.16 x 10 5 

192 

3 

4 

48 

1.60 x 10 5 

192 

7 

8 

24 

3.34 x 10 5 

192 

15 

16 

12 

unstable 

256 

1 

2 

128 

5.28 x 10 5 

256 

3 

4 

64 

7.27 x 10 5 

256 

7 

8 

32 

1.53 x 10 6 

256 

15 

16 

16 

3.06 x 10 6 

384 

1 

2 

192 

2.84 x 10 6 

384 

3 

4 

96 

3.87 x 10 6 

384 

7 

8 

48 

8.11 x 10 6 

384 

15 

16 

24 

1.16 x 10 7 * 


(*Case stopped at t = 10.75) 


Table 3: Case summary for the Taylor-Green vortex evolution, Mq = 0.1, Re = 1600. 


Results 

For each run the temporal evolution of the kinetic energy 

E k = 71 / \p w ' ( 8 ) 

“ Jn 

was monitored. The evolution of the kinetic energy dissipation rate e = was computed using a finite difference 

approximation. Figure 3 plots the dimensionless kinetic energy dissipation rate e vs time for 2nd- ( p = 1), 4th- 
(p = 3), 8th- (p = 7), 16th- (p = 15) order schemes. Figure 3 also plots the dissipation computed for an incompressible 
simulation using a spectral code on a 512 3 grid [7]. For the 2nd-order scheme there is significant variation in the 
dissipation rate between the four sets of meshes, as well as the reference incompressible spectral data. For 4th-order 
solutions, the computed results are converging to the spectral data and there is little noticeable difference between 
the three finest mesh sizes on the given scale. At 8th and 16th order, there is insufficient numerical dissipation on 
the coarsest mesh leading to instability. Using 192 degrees of freedom in each coordinate direction, the 8th-order 
scheme is stable while the 16th-order scheme becomes unstable shortly after the point of peak dissipation. We note, 
however, that up to the point of failure there is still very good matching with the spectral data. A close-up view of 
the evolution of the dissipation rate at the point of peak dissipation is given in Figure 4 for 4th-, 8th- and 16th-order 
solutions. 

We assess the quality of our numerical solutions by computing individual terms in the kinetic energy evolution 
equation. For incompressible flow the kinetic energy dissipation rate is equal to 2/i£, where £ is the enstrophy, 
computed as: 


£ 


1 

n 


• l od£l 


(9) 


where uj = V x v is the vorticity. 


For compressible flow, the kinetic energy dissipation rate is given by the sum of 
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(a) 2nd order (p = 1) (b) 4th order (p = 3) 




(c) 8th order (p = 7) (d) 16th order (p = 15) 

Figure 3: Kinetic energy dissipation rate for the Taylor-Green vortex evolution, M 0 = 0.1, Re = 1600. 
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(a) 4th order ( p = 3) (b) 8th order ( p = 7) 



(c) 16th order (p = 15) 


Figure 4: Kinetic energy dissipation rate at peak dissipation for the Taylor-Green vortex evolution, M 0 = 0.1, 
Re = 1600. 
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three contributions e = e\ + e 2 + £3 = — : 


£i = 

j- [ 2/i S : S dtt 

(10) 

£2 = 

j- 1 A (V • v) 2 dfl 
J n 

(11) 

£3 = 

— — / p v • v dn 
- ' J n 

(12) 


where S = |(Vv + Vv T ) is the strain rate tensor. We note that £, ei, e 2 and e 3 are computed using the “lifted” 
gradients in order to be consistent with our DG discretization. 

Since the flow is nearly incompressible, we expect that the dissipation due to the bulk viscosity (e 2 ) and the pressure 
strain term (e 3 ) to be small. The kinetic energy dissipation rate is then approximately equal to e k 2pE ss e±. 
Differences between these quantities indicates the presence of compressibility effects and numerical dissipation. 



Time 


(a) 2nd order ( p = 1), 192 3 



Time 


(b) 2nd order (p = 1), 384 3 


Figure 5: Evolution of terms in kinetic energy balance equation for the Taylor-Green vortex evolution, Mo = 0.1, 
Re = 1600. 


Figure 5 plots the temporal evolution of e, ei, e 2 and e 3 for the 2nd-order (p = 1) simulations. It appears as 
though there is a significant contribution to the dissipation due to the pressure strain term, e 3 . A significant amount 
of numerical dissipation is present as indicated by the difference between e and £1 + £2 + £ 3 - Figure 6 shows the 
corresponding plots for the 4th ( p = 3) and 8th {p = 4) simulations with 192 3 degrees of freedom. For these 
simulations there is much less contribution of the pressure strain term and less numerical dissipation. Figure 7 shows 
the pressure strain term through the sequence of mesh refinements. 

Finally we plot contours of vorticity magnitude on the face x = —nL at t = 8. Figure 8 plots the vorticity magnitude 
on the finest mesh for each polynomial order. 
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Figure 6: Evolution of terms in kinetic energy balance equation for the Taylor-Green vortex evolution, Mq = 0.1, 
Re = 1600. 
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(d) 16th order (p = 15) 


Figure 7: Pressure- strain dissipation for the Taylor- Green vortex evolution, Mq = 0.1, Re = 1600. 
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(a) 2nd order ( p = 1) (b) 4th order ( p = 3) 



(c) 8th order (p = 7) (d) 16th order ( p = 15) 


Figure 8: Contours of vorticity magnitude at x = — nL , t = 0 for the Taylor-Green vortex evolution, M 0 = 0.1, 
Re = 1600 using 256 x 256 x 256 degrees of freedom. 


10 






Case summaries for 2nd International Workshop on Higher-Order CFD Methods, May 27-28, 2013, Koln, Germany 


[3] Gassner, G. and Kopriva, D. A., “A comparison of the dispersion and dissipation errors of Gauss and Gauss- 
Lobatto discontinuous Galerkin spectral element methods,” SIAM J. Sci. Comput ., Vol. 33, 2011, pp. 2560-2579. 

[4] Roe, P. L., “Approximate Riemann solvers, parameter vectors, and difference schemes,” Journal of Computational 
Physics, Vol. 43, No. 2, 1981, pp. 357 372. 

[5] Bassi, F. and Rebay, S., “GMRES discontinuous Galerkin solution of the compressible Navier-Stokes equations,” 
Discontinuous Galerkin Methods: Theory, Computation and Applications, edited by K. Cockburn and Shu, 
Springer, Berlin, 2000, pp. 197-208. 

[6] Wang, Z. J., “1st International Workshop on High-Order CFD Methods,” 
http://zjw.public.iastate.edu/hiocfd.html, 2012. 

[7] van Ress, W., Leonard, A., Pullin, D., and Koumoutsakos, P., “A comparison of vortex and pseudo-spectral 
methods for the simulation of periodic vortical flows at high Reynolds number,” Journal of Computational 
Physics, Vol. 230, 2011, pp. 2794-2805. 


11 



